
/* 
* sepacf.c - This file is a part of the sep-library 
*
* Copyright (C) 2008 Jesper Schmidt Hansen 
* 
* License: GPL - see COPYING for copying conditions.
* There is ABSOLUTELY NO WARRANTY, not even for MERCHANTIBILITY or
* FITNESS FOR A PARTICULAR PURPOSE.
*
* Contact: schmidt@zigzak.net
*/

#include "sepacf.h"

char *sep_acaf_warning( void ){
  char *retval;
  
  retval = "All elements in vacf array not called the same number of times";
  
  return retval;
}

// Atomic vel. autocorrelations
double sep_eval_avacf(seppart *ptr, char type, double **vel0, 
                     int opt, int npart, int ndim){
  double retval = 0.0;
  int n, k;

  switch ( opt ){
    case 0:
      for ( n=0; n<npart; n++ ){
        if ( ptr[n].type == type ){
          for ( k=0; k<ndim; k++ )
            vel0[n][k] = ptr[n].v[k];
        }
      }
      break;
      
    case 1:
      for ( n=0; n<npart; n++ ){
        if ( ptr[n].type == type ) {
          for ( k=0; k<ndim; k++ )
            retval += ptr[n].v[k]*vel0[n][k]; 
        }
      }         
      retval /= npart*ndim;
      break;
  }

  return retval;
}


int sep_sample_avacf(seppart *ptr, char type, double *vac, 
                     unsigned int lvec, long unsigned int *ncall,
                     int opt, int npart, int ndim){
  static unsigned count = 0;
  static double **vel_0;
  int ntype = 0, retval = 0;
  unsigned int n;
  
  switch ( opt ) {
    case 0:
    
    ntype = sep_count_type(ptr, type, npart);               
    if ( count == 0 ){
      vel_0 = sep_matrix(ntype, ndim);
      sep_eval_avacf(ptr, type, vel_0, 0, npart, ndim);
    }
    vac[count] += sep_eval_avacf(ptr, type, vel_0, 1, npart, ndim);

    count++;
    if ( count == lvec ){
      sep_free_matrix(vel_0, ntype);
      count = 0;
      (*ncall)++;
      retval = 1;
    }
   
    break;
	 
  case 1:
    
    if ( count != 0 ){
      sep_free_matrix(vel_0, ntype);
      sep_warning("%s at line %d: %s \n",
                  sep_acaf_warning(), __func__, __LINE__);
    }
			
    for ( n=0; n<lvec; n++ ){
      vac[n] /= (*ncall);
    }
			
    break;
  }

  return retval;
}


// Atomic press. tensor autocorr
double sep_apacf_kin(seppart *ptr, char type, int npart){
  double retval=0.0;
  int n;

  for ( n=0; n<npart; n++ ){
    if ( ptr[n].type == type )
      retval += ptr[n].v[0]*ptr[n].v[1];
  }
  
  return retval;

}


double sep_2D_apacf_nonbon(seppart *ptr, char type, double (*fun)(double, char),
                          double cf, sep2D sys){
  int n, i1, i2, k;
  double retval, r[3], r2, cf2 = cf*cf, ft;

  retval = 0.0;
  for ( i1=0; i1<sys.npart; i1++ ){
    if ( ptr[i1].type != type ) continue;
    n = 0;
    while ( 1 ){
      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 || ptr[i2].type != type ) break; 
      r2 = 0.0;
      for ( k=0; k<2; k++ ){
	r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	if ( sys.bound[k] == 'p' )
	  sep_Wrap(r[k], sys.length[k]);
	r2 += r[k]*r[k];
      }
      if ( r2 < cf2 ){
	ft = fun(r2, 'f');
	retval += r[0]*ft*r[1];
      }
      n++;
    }
  }
  
  return retval;
}


double sep_3D_apacf_nonbon(seppart *ptr, char type, double (*fun)(double, char),
                           double cf, sep3D sys){ 
  int n, i1, i2, k;
  double retval, r[3], r2, cf2 = cf*cf, ft;

  retval = 0.0;
  for ( i1=0; i1<sys.npart; i1++ ){
    if ( ptr[i1].type != type ) continue;
    n = 0;
    while ( 1 ){
      i2 = ptr[i1].neighb[n];
      if ( i2 == -1 || ptr[i2].type != type ) break; 
      r2 = 0.0;
      for ( k=0; k<3; k++ ){
	r[k] = ptr[i1].x[k] - ptr[i2].x[k];
	if ( sys.bound[k] == 'p' )
	  sep_Wrap(r[k], sys.length[k]);
	r2 += r[k]*r[k];
      }
      if ( r2 < cf2 ){
	ft = fun(r2, 'f');
	retval += r[0]*ft*r[1];
      }
      n++;
    }
  }
  
  return retval;
}


void sep_3D_eval_apress(seppart *ptr, double **P, char type, 
			double (*fun)(double, char),
			double cf, sep3D sys){
  long int i, j, k, kk, n;
  double pekin[3][3], ppot[3][3], r2, r[3], ft;
  const double cf2 = cf*cf;
  
  for ( k=0; k<3; k++ )
    for ( kk=0; kk<3; kk++ )
      pekin[k][kk] = ppot[k][kk] = 0.0;
  
  // Kinetic part
  for ( i=0; i<sys.npart; i++ ){
    if ( ptr[i].type == type ){
      for ( k=0; k<3; k++ )
	for ( kk=0; kk<3; kk++ )
	  pekin[k][kk] += ptr[i].m*ptr[i].v[k]*ptr[i].v[kk];
    }
  }	
  
  // Potential part
  for ( i=0; i<sys.npart; i++ ){
    if ( ptr[i].type != type ) continue;
    n = 0;
    while ( 1 ){
      j = ptr[i].neighb[n];
      if ( j == -1 || ptr[j].type != type ) break; 
      r2 = 0.0;
      for ( k=0; k<3; k++ ){
	r[k] = ptr[i].x[k] - ptr[j].x[k];
	if ( sys.bound[k] == 'p' )
	  sep_Wrap(r[k], sys.length[k]);
	r2 += r[k]*r[k];
      }
      if ( r2 < cf2 ){
	ft = fun(r2, 'f');
	for ( k=0; k<3; k++ )
	  for ( kk=0; kk<3; kk++ )
	    ppot[k][kk] += r[k]*ft*r[kk];
      }
      n++;	
    }
  }
  
  // Adding the two contributions together
  for ( k=0; k<3; k++ )
    for ( kk=0; kk<3; kk++ )
      P[k][kk] = (pekin[k][kk] + ppot[k][kk])/sys.volume; 	
  
}



void sep_3D_eval_mpress(seppart *ptr, double **P, char type, 
                        double (*fun)(double, char), const double cf,  
                        int nbeads, int npoly, sep3D sys){
  double **cm, **pkin, **ppot, veli[3], rij[3], fij[3], rinjm[3], rsq, ft;
  int i, j, n, m, in, jm, k, kk, *heads;
  const double cfsq = cf*cf;
  
  // Initialising arrays
  pkin = sep_matrix(3,3);
  ppot = sep_matrix(3,3);
  cm = sep_matrix(npoly, 3);  
  heads = sep_vector_int(npoly);
  
  // Kinetic contribution
  sep_linmol_heads(ptr, heads, type, nbeads, sys.npart);
  for ( i=0; i<npoly; i++ ){

    sep_vector_set(veli, 3, 0.0);
    for ( n=0; n<nbeads; n++ ){
      in = heads[i] + n;
      for ( k=0; k<3; k++ )
        veli[k] += ptr[in].v[k];
    }
    
    for ( k=0; k<3; k++ )    
      for ( kk=0; kk<3; kk++ )
        pkin[k][kk] += veli[k]*veli[kk]/nbeads;
  }
  
  // Potential part - brute force  
  sep_3D_linmol_cm(ptr, cm, type, nbeads, sys);
  for ( i=0; i<npoly; i++ ){
    for ( j=0; j<npoly; j++ ){
      
      if ( i==j ) continue;
      
      for ( k=0; k<3; k++ ){
        rij[k] = cm[i][k] - cm[j][k];
        sep_Wrap( rij[k], sys.length[k] );  
      }
      
      sep_vector_set(fij, 3, 0.0);   
      for ( n=0; n<nbeads; n++) {
          in = heads[i] + n;	
        for ( m=0; m<nbeads; m++ ) {
          jm = heads[j] + m;	
          rsq = 0.0;
          for ( k=0; k<3; k++ ){
            rinjm[k] = ptr[in].x[k]-ptr[jm].x[k];
            sep_Wrap(rinjm[k], sys.length[k]);
            rsq += rinjm[k]*rinjm[k];
          }
          if ( rsq < cfsq ){
            ft = fun(rsq, 'f');  
            for ( k=0; k<3; k++ )
              fij[k] += ft*rinjm[k];
          }
        }
      }
      
      for ( k=0; k<3; k++ )
        for ( kk=0; kk<3; kk++ )
          ppot[k][kk] += rij[k]*fij[kk];
      
    }
  } 
  
  // Adding the kinetic and potential contriubtions
  for ( k=0; k<3; k++ ) {
    for ( kk=0; kk<3; kk++ ) {
      P[k][kk] = pkin[k][kk] + 0.5*ppot[k][kk];
      //P[k][kk] = pkin[k][kk];
      //P[k][kk] = 0.5*ppot[k][kk];
    }
  }
      
  free(heads);             
  sep_free_matrix(ppot, 3);
  sep_free_matrix(pkin, 3);   
  sep_free_matrix(cm, npoly);
   
}


int sep_3D_sample_apacf(seppart *ptr, char type, double *pxy, int lvec,
                         double (*fun)(double, char), double cf,
                         unsigned long int *ncall, int opt, sep3D sys) { 
  int n, m, retval = 0;
  static int count = 0;
  static double *pxy_0;
    
  switch (opt) {
  case 0:

    if ( count == 0 ) 
      pxy_0 = sep_vector(lvec);
       
    pxy_0[count] = sep_apacf_kin(ptr, type, sys.npart);
    pxy_0[count] += sep_3D_apacf_nonbon(ptr, type, fun, cf, sys);
    count++;
 
    if ( count == lvec ){
      for ( n=0; n<lvec; n++ ){
	for ( m=0; m<lvec-n; m++ ){
	  pxy[n] += pxy_0[m]*pxy_0[m+n];
	}
      }
      free(pxy_0);
      count = 0;
      (*ncall)++;
      retval = 1;
    }
    
    break;
  case 1:
    
    if ( count != 0 )
      sep_warning("%s at line %d: %s\n", 
                  sep_acaf_warning, __func__, __LINE__);
      
    for ( n=0; n<lvec; n++ ){
      pxy[n] /= (*ncall)*(lvec-n);
    }

    break;
  }

  return retval;  
}



int sep_2D_sample_apacf(seppart *ptr, char type, double *pxy, int lvec,
                         double (*fun)(double, char), double cf,
                         unsigned long int *ncall, int opt, sep2D sys) { 
  int n, m, retval = 0;
  static int count = 0;
  static double *pxy_0;
    
  switch (opt) {
  case 0:

    if ( count == 0 ) 
      pxy_0 = sep_vector(lvec);
       
    pxy_0[count] = sep_apacf_kin(ptr, type, sys.npart);
    pxy_0[count] += sep_2D_apacf_nonbon(ptr, type, fun, cf, sys);
    count++;
 
    if ( count == lvec ){
      for ( n=0; n<lvec; n++ ){
	for ( m=0; m<lvec-n; m++ ){
	  pxy[n] += pxy_0[m]*pxy_0[m+n];
	}
      }
      free(pxy_0);
      count = 0;
      (*ncall)++;
      retval = 1;
    }
    
    break;
  case 1:
    
    if ( count != 0 )
      sep_warning("%s at line %d: %s\n", sep_acaf_warning(), 
                  __func__, __LINE__);
      
    for ( n=0; n<lvec; n++ ){
      pxy[n] /= (*ncall)*(lvec-n);
    }

    break;
  }

  return retval;  
}
 
double sep_3D_sample_mpacf(seppart *ptr, char type, 
                        double *pacf_os, double *pacf_as, int lvec,  
                        double (*fun)(double, char), const double cf,
                        int npoly, int nbeads, unsigned long int *ncall, 
                        int opt, sep3D sys){ 
  int n, m;
  static int count = 0, count1=0;
  static double ***mol_press_asym, ***mol_press_sym,
                **accum_sym, **accum_asym;
  double **P, trace, press=0.0;
  const double vol = sys.length[0]*sys.length[1]*sys.length[2];
  
  if ( count1 == 0 ) {
    count1 = 1;
    accum_sym = sep_matrix(lvec, 9);
    accum_asym = sep_matrix(lvec, 3);
  }
  
  switch (opt) {
  case 0:
       
    if ( count == 0 ){
      mol_press_asym = sep_tensor(3, 3, lvec);
      mol_press_sym  = sep_tensor(3, 3, lvec);
    } 
      
    P = sep_matrix(3, 3);
    sep_3D_eval_mpress(ptr, P, type, fun, cf, nbeads, npoly, sys);
    trace = sep_trace(P, 3);
    press = trace/(vol*3.0);
   
    for ( n=0; n<3; n++ ){
      for ( m=0; m<3; m++ ){
	mol_press_sym[n][m][count] = 0.5*(P[n][m] + P[m][n]);
        if ( n==m )
          mol_press_sym[n][m][count] -= trace/3.0;
	mol_press_asym[n][m][count] = 0.5*(P[n][m] - P[m][n]);     
      }
    }     
    sep_free_matrix(P,3);       
        
    count++;
    if ( count == lvec ){
      for ( n=0; n<lvec; n++ ){
        for ( m=0; m<lvec-n; m++ ){
	  accum_sym[n][0] +=mol_press_sym[0][0][m]*mol_press_sym[0][0][m+n];
          accum_sym[n][1] +=mol_press_sym[0][1][m]*mol_press_sym[0][1][m+n];
          accum_sym[n][2] +=mol_press_sym[0][2][m]*mol_press_sym[0][2][m+n];
          accum_sym[n][3] +=mol_press_sym[1][1][m]*mol_press_sym[1][1][m+n];
          accum_sym[n][4] +=mol_press_sym[1][2][m]*mol_press_sym[1][2][m+n];
          accum_sym[n][5] +=mol_press_sym[2][2][m]*mol_press_sym[2][2][m+n];
          accum_sym[n][6] +=mol_press_sym[1][0][m]*mol_press_sym[1][0][m+n];
	  accum_sym[n][7] +=mol_press_sym[2][0][m]*mol_press_sym[2][0][m+n];
          accum_sym[n][8] +=mol_press_sym[2][1][m]*mol_press_sym[2][1][m+n];
         
          accum_asym[n][0]+=mol_press_asym[0][1][m]*mol_press_asym[0][1][m+n];
          accum_asym[n][1]+=mol_press_asym[2][0][m]*mol_press_asym[2][0][m+n];
          accum_asym[n][2]+=mol_press_asym[1][2][m]*mol_press_asym[1][2][m+n];
        }
      }
      count = 0;
      (*ncall)++;
      sep_free_tensor(mol_press_sym, 3, 3);
      sep_free_tensor(mol_press_asym, 3, 3);
    }
 
    break;
  case 1:
    
    if ( count != 0 ){
       sep_warning("%s at line %d: %s\n", sep_acaf_warning(), 
                    __func__, __LINE__);
    }
    
    for ( n=0; n<lvec; n++ ){
      pacf_os[n] = 0.0;
      pacf_as[n] = 0.0;	
      for ( m=0; m<9; m++ ){
        pacf_os[n] += accum_sym[n][m];
        if ( m<3 )
          pacf_as[n] += accum_asym[n][m];
      }
      pacf_os[n] /= (*ncall)*10.0*(lvec-n);
      pacf_as[n] /= (*ncall)*3.0*(lvec-n);
    }
    
    sep_free_matrix(accum_sym, lvec);
    sep_free_matrix(accum_asym, lvec);
    
    break;
  }

  return press;  
}

